"""PPG处理结果分析和可视化

提供PPG处理结果的详细分析和可视化功能，包括：
1. 结果统计分析
2. 信号和峰值可视化
3. HRV分析图表
4. 质量评估报告
5. 结果对比分析

作者: PPG算法包开发团队
版本: 2.0.0
"""

import numpy as np
import logging
from typing import List, Dict, Any, Optional, Tuple
import json
from pathlib import Path

try:
    from ..core.data_types import PPGResults, HRVResults, QualityResults
except ImportError:
    # 处理直接运行时的导入问题
    import sys
    sys.path.append(str(Path(__file__).parent.parent))
    from core.data_types import PPGResults, HRVResults, QualityResults

# 设置日志
logging.basicConfig(level=logging.INFO)
logger = logging.getLogger(__name__)

# 尝试导入可视化库
try:
    import matplotlib.pyplot as plt
    import matplotlib.patches as patches
    from matplotlib.gridspec import GridSpec
    MATPLOTLIB_AVAILABLE = True
except ImportError:
    MATPLOTLIB_AVAILABLE = False
    logger.warning("matplotlib未安装，可视化功能将不可用")

try:
    import seaborn as sns
    SEABORN_AVAILABLE = True
    if MATPLOTLIB_AVAILABLE:
        sns.set_style("whitegrid")
except ImportError:
    SEABORN_AVAILABLE = False


def analyze_ppg_results(results: PPGResults, detailed: bool = True) -> Dict[str, Any]:
    """分析PPG处理结果
    
    参数:
        results: PPG处理结果
        detailed: 是否进行详细分析
        
    返回:
        Dict: 分析结果字典
    """
    print("=" * 60)
    print("PPG处理结果分析")
    print("=" * 60)
    
    analysis = {
        'basic_info': analyze_basic_info(results),
        'signal_analysis': analyze_signal_characteristics(results),
        'peak_analysis': analyze_peak_characteristics(results),
        'heart_rate_analysis': analyze_heart_rate(results),
        'quality_analysis': analyze_quality(results)
    }
    
    if results.hrv_results and detailed:
        analysis['hrv_analysis'] = analyze_hrv_detailed(results.hrv_results)
    
    # 打印分析结果
    print_analysis_summary(analysis)
    
    return analysis


def analyze_basic_info(results: PPGResults) -> Dict[str, Any]:
    """分析基本信息
    
    参数:
        results: PPG处理结果
        
    返回:
        Dict: 基本信息分析
    """
    sampling_rate = len(results.original_signal) / (len(results.original_signal) / 100.0)  # 估算采样率
    duration = len(results.original_signal) / sampling_rate
    
    return {
        'pipeline_name': results.pipeline_name,
        'signal_length': len(results.original_signal),
        'processed_length': len(results.processed_signal),
        'estimated_sampling_rate': sampling_rate,
        'estimated_duration': duration,
        'peaks_detected': len(results.peaks),
        'processing_steps': len(results.processing_info) if results.processing_info else 0
    }


def analyze_signal_characteristics(results: PPGResults) -> Dict[str, Any]:
    """分析信号特征
    
    参数:
        results: PPG处理结果
        
    返回:
        Dict: 信号特征分析
    """
    original = results.original_signal
    processed = results.processed_signal
    
    # 原始信号统计
    orig_stats = {
        'mean': float(np.mean(original)),
        'std': float(np.std(original)),
        'min': float(np.min(original)),
        'max': float(np.max(original)),
        'range': float(np.ptp(original)),
        'skewness': float(calculate_skewness(original)),
        'kurtosis': float(calculate_kurtosis(original))
    }
    
    # 处理后信号统计
    proc_stats = {
        'mean': float(np.mean(processed)),
        'std': float(np.std(processed)),
        'min': float(np.min(processed)),
        'max': float(np.max(processed)),
        'range': float(np.ptp(processed)),
        'skewness': float(calculate_skewness(processed)),
        'kurtosis': float(calculate_kurtosis(processed))
    }
    
    # 信号改善指标
    snr_improvement = calculate_snr_improvement(original, processed)
    
    return {
        'original_signal': orig_stats,
        'processed_signal': proc_stats,
        'snr_improvement': snr_improvement,
        'processing_effect': {
            'mean_change': proc_stats['mean'] - orig_stats['mean'],
            'std_change': proc_stats['std'] - orig_stats['std'],
            'range_reduction': (orig_stats['range'] - proc_stats['range']) / orig_stats['range']
        }
    }


def analyze_peak_characteristics(results: PPGResults) -> Dict[str, Any]:
    """分析峰值特征
    
    参数:
        results: PPG处理结果
        
    返回:
        Dict: 峰值特征分析
    """
    if len(results.peaks) == 0:
        return {'error': 'no_peaks_detected'}
    
    peaks = results.peaks
    signal = results.processed_signal
    
    # 峰值间隔分析
    intervals = np.diff(peaks)
    
    # 峰值高度分析
    peak_heights = signal[peaks]
    
    # 峰值规律性分析
    regularity = calculate_peak_regularity(intervals)
    
    return {
        'peak_count': len(peaks),
        'intervals': {
            'mean': float(np.mean(intervals)),
            'std': float(np.std(intervals)),
            'min': float(np.min(intervals)),
            'max': float(np.max(intervals)),
            'cv': float(np.std(intervals) / np.mean(intervals)) if np.mean(intervals) > 0 else 0
        },
        'heights': {
            'mean': float(np.mean(peak_heights)),
            'std': float(np.std(peak_heights)),
            'min': float(np.min(peak_heights)),
            'max': float(np.max(peak_heights))
        },
        'regularity_score': regularity,
        'peak_density': len(peaks) / len(signal)  # 峰值密度
    }


def analyze_heart_rate(results: PPGResults) -> Dict[str, Any]:
    """分析心率特征
    
    参数:
        results: PPG处理结果
        
    返回:
        Dict: 心率分析
    """
    if results.heart_rates is None or len(results.heart_rates) == 0:
        return {'error': 'no_heart_rate_data'}
    
    hr = results.heart_rates
    
    # 心率变异性指标
    hr_variability = {
        'range': float(np.max(hr) - np.min(hr)),
        'std': float(np.std(hr)),
        'cv': float(np.std(hr) / np.mean(hr)) if np.mean(hr) > 0 else 0
    }
    
    # 心率趋势分析
    trend = analyze_heart_rate_trend(hr)
    
    # 异常心率检测
    outliers = detect_heart_rate_outliers(hr)
    
    return {
        'mean_heart_rate': results.mean_heart_rate,
        'heart_rate_stats': {
            'min': float(np.min(hr)),
            'max': float(np.max(hr)),
            'median': float(np.median(hr)),
            'q25': float(np.percentile(hr, 25)),
            'q75': float(np.percentile(hr, 75))
        },
        'variability': hr_variability,
        'trend': trend,
        'outliers': outliers,
        'stability': classify_heart_rate_stability(hr_variability['cv'])
    }


def analyze_quality(results: PPGResults) -> Dict[str, Any]:
    """分析信号质量
    
    参数:
        results: PPG处理结果
        
    返回:
        Dict: 质量分析
    """
    if not results.quality_results:
        return {'error': 'no_quality_data'}
    
    quality = results.quality_results
    
    # 质量等级解释
    grade_explanation = {
        'excellent': '信号质量优秀，适合所有类型的分析',
        'good': '信号质量良好，适合大多数分析',
        'fair': '信号质量一般，建议谨慎使用分析结果',
        'poor': '信号质量较差，不建议用于精确分析'
    }
    
    return {
        'overall_score': quality.quality_score,
        'quality_grade': quality.quality_grade,
        'grade_explanation': grade_explanation.get(quality.quality_grade, '未知等级'),
        'snr': quality.snr,
        'peak_regularity': quality.peak_regularity,
        'artifacts_ratio': quality.artifacts_ratio,
        'quality_recommendations': generate_quality_recommendations(quality)
    }


def analyze_hrv_detailed(hrv: HRVResults) -> Dict[str, Any]:
    """详细分析HRV结果
    
    参数:
        hrv: HRV分析结果
        
    返回:
        Dict: 详细HRV分析
    """
    # HRV指标解释
    interpretations = {
        'rmssd': interpret_rmssd(hrv.rmssd),
        'sdnn': interpret_sdnn(hrv.sdnn),
        'pnn50': interpret_pnn50(hrv.pnn50),
        'lf_hf_ratio': interpret_lf_hf_ratio(hrv.lf_hf_ratio)
    }
    
    # 自主神经系统状态评估
    ans_assessment = assess_autonomic_nervous_system(hrv)
    
    return {
        'time_domain': {
            'rmssd': hrv.rmssd,
            'sdnn': hrv.sdnn,
            'pnn50': hrv.pnn50,
            'mean_rr': hrv.mean_rr,
            'std_rr': hrv.std_rr
        },
        'frequency_domain': {
            'lf_power': hrv.lf_power,
            'hf_power': hrv.hf_power,
            'lf_hf_ratio': hrv.lf_hf_ratio,
            'total_power': hrv.total_power
        },
        'nonlinear': {
            'sd1': hrv.sd1,
            'sd2': hrv.sd2,
            'sd1_sd2_ratio': hrv.sd1 / hrv.sd2 if hrv.sd2 > 0 else 0
        },
        'interpretations': interpretations,
        'ans_assessment': ans_assessment,
        'rr_intervals_count': len(hrv.rr_intervals)
    }


def visualize_ppg_results(results: PPGResults, save_path: Optional[str] = None, 
                         show_plot: bool = True) -> bool:
    """可视化PPG处理结果
    
    参数:
        results: PPG处理结果
        save_path: 保存路径（可选）
        show_plot: 是否显示图表
        
    返回:
        bool: 是否成功创建可视化
    """
    if not MATPLOTLIB_AVAILABLE:
        print("matplotlib未安装，无法创建可视化")
        return False
    
    try:
        # 创建子图布局
        fig = plt.figure(figsize=(15, 12))
        gs = GridSpec(4, 2, figure=fig, hspace=0.3, wspace=0.3)
        
        # 1. 原始信号和处理后信号对比
        ax1 = fig.add_subplot(gs[0, :])
        plot_signal_comparison(ax1, results)
        
        # 2. 峰值检测结果
        ax2 = fig.add_subplot(gs[1, :])
        plot_peak_detection(ax2, results)
        
        # 3. 心率变化
        ax3 = fig.add_subplot(gs[2, 0])
        plot_heart_rate(ax3, results)
        
        # 4. HRV分析（如果有）
        ax4 = fig.add_subplot(gs[2, 1])
        if results.hrv_results:
            plot_hrv_poincare(ax4, results.hrv_results)
        else:
            ax4.text(0.5, 0.5, 'HRV数据不可用', ha='center', va='center', transform=ax4.transAxes)
            ax4.set_title('HRV Poincaré图')
        
        # 5. 质量评估
        ax5 = fig.add_subplot(gs[3, 0])
        plot_quality_assessment(ax5, results)
        
        # 6. 频谱分析
        ax6 = fig.add_subplot(gs[3, 1])
        plot_frequency_spectrum(ax6, results)
        
        plt.suptitle(f'PPG处理结果分析 - {results.pipeline_name}', fontsize=16, fontweight='bold')
        
        if save_path:
            plt.savefig(save_path, dpi=300, bbox_inches='tight')
            print(f"可视化结果已保存到: {save_path}")
        
        if show_plot:
            plt.show()
        else:
            plt.close()
        
        return True
        
    except Exception as e:
        logger.error(f"可视化创建失败: {e}")
        return False


def plot_signal_comparison(ax, results: PPGResults):
    """绘制信号对比图"""
    t_orig = np.arange(len(results.original_signal))
    t_proc = np.arange(len(results.processed_signal))
    
    ax.plot(t_orig, results.original_signal, 'lightgray', alpha=0.7, label='原始信号')
    ax.plot(t_proc, results.processed_signal, 'blue', linewidth=1.5, label='处理后信号')
    
    ax.set_title('信号预处理对比')
    ax.set_xlabel('采样点')
    ax.set_ylabel('幅值')
    ax.legend()
    ax.grid(True, alpha=0.3)


def plot_peak_detection(ax, results: PPGResults):
    """绘制峰值检测结果"""
    t = np.arange(len(results.processed_signal))
    ax.plot(t, results.processed_signal, 'blue', linewidth=1, label='处理后信号')
    
    if len(results.peaks) > 0:
        ax.plot(results.peaks, results.processed_signal[results.peaks], 
                'ro', markersize=6, label=f'检测峰值 ({len(results.peaks)}个)')
    
    ax.set_title('峰值检测结果')
    ax.set_xlabel('采样点')
    ax.set_ylabel('幅值')
    ax.legend()
    ax.grid(True, alpha=0.3)


def plot_heart_rate(ax, results: PPGResults):
    """绘制心率变化"""
    if results.heart_rates is not None and len(results.heart_rates) > 0:
        t = np.arange(len(results.heart_rates))
        ax.plot(t, results.heart_rates, 'g-', linewidth=2, marker='o', markersize=4)
        ax.axhline(y=results.mean_heart_rate, color='r', linestyle='--', 
                  label=f'平均心率: {results.mean_heart_rate:.1f} BPM')
        ax.fill_between(t, results.heart_rates, alpha=0.3, color='green')
    else:
        ax.text(0.5, 0.5, '心率数据不可用', ha='center', va='center', transform=ax.transAxes)
    
    ax.set_title('心率变化')
    ax.set_xlabel('心跳序号')
    ax.set_ylabel('心率 (BPM)')
    ax.legend()
    ax.grid(True, alpha=0.3)


def plot_hrv_poincare(ax, hrv: HRVResults):
    """绘制HRV Poincaré图"""
    if len(hrv.rr_intervals) < 2:
        ax.text(0.5, 0.5, 'RR间期数据不足', ha='center', va='center', transform=ax.transAxes)
        return
    
    rr1 = hrv.rr_intervals[:-1]
    rr2 = hrv.rr_intervals[1:]
    
    ax.scatter(rr1, rr2, alpha=0.6, s=20)
    
    # 添加SD1和SD2椭圆
    if hrv.sd1 > 0 and hrv.sd2 > 0:
        center = (np.mean(rr1), np.mean(rr2))
        ellipse = patches.Ellipse(center, 2*hrv.sd2, 2*hrv.sd1, 
                                angle=45, fill=False, color='red', linewidth=2)
        ax.add_patch(ellipse)
    
    ax.set_title(f'HRV Poincaré图\nSD1={hrv.sd1:.1f}, SD2={hrv.sd2:.1f}')
    ax.set_xlabel('RR(n) (ms)')
    ax.set_ylabel('RR(n+1) (ms)')
    ax.grid(True, alpha=0.3)


def plot_quality_assessment(ax, results: PPGResults):
    """绘制质量评估"""
    if not results.quality_results:
        ax.text(0.5, 0.5, '质量数据不可用', ha='center', va='center', transform=ax.transAxes)
        return
    
    quality = results.quality_results
    
    # 创建质量指标条形图
    metrics = ['总体评分', '信噪比', '峰值规律性']
    values = [quality.quality_score, quality.snr/20, quality.peak_regularity]  # 归一化SNR
    colors = ['blue', 'green', 'orange']
    
    bars = ax.bar(metrics, values, color=colors, alpha=0.7)
    
    # 添加数值标签
    for bar, value in zip(bars, [quality.quality_score, quality.snr, quality.peak_regularity]):
        height = bar.get_height()
        ax.text(bar.get_x() + bar.get_width()/2., height + 0.01,
                f'{value:.2f}', ha='center', va='bottom')
    
    ax.set_title(f'信号质量评估 - {quality.quality_grade}')
    ax.set_ylabel('评分')
    ax.set_ylim(0, 1.2)
    ax.grid(True, alpha=0.3, axis='y')


def plot_frequency_spectrum(ax, results: PPGResults):
    """绘制频谱分析"""
    try:
        from scipy.signal import welch
        
        # 估算采样率
        sampling_rate = 100.0  # 默认采样率
        
        freqs, psd = welch(results.processed_signal, fs=sampling_rate, nperseg=min(256, len(results.processed_signal)//4))
        
        ax.semilogy(freqs, psd)
        ax.set_title('功率谱密度')
        ax.set_xlabel('频率 (Hz)')
        ax.set_ylabel('功率谱密度')
        ax.grid(True, alpha=0.3)
        
        # 标记心率频率范围
        hr_freq_min = 0.5  # 30 BPM
        hr_freq_max = 3.0  # 180 BPM
        ax.axvspan(hr_freq_min, hr_freq_max, alpha=0.2, color='red', label='心率频率范围')
        ax.legend()
        
    except ImportError:
        ax.text(0.5, 0.5, 'scipy不可用', ha='center', va='center', transform=ax.transAxes)


# 辅助函数
def calculate_skewness(data: np.ndarray) -> float:
    """计算偏度"""
    mean = np.mean(data)
    std = np.std(data)
    if std == 0:
        return 0
    return np.mean(((data - mean) / std) ** 3)


def calculate_kurtosis(data: np.ndarray) -> float:
    """计算峰度"""
    mean = np.mean(data)
    std = np.std(data)
    if std == 0:
        return 0
    return np.mean(((data - mean) / std) ** 4) - 3


def calculate_snr_improvement(original: np.ndarray, processed: np.ndarray) -> float:
    """计算信噪比改善"""
    # 简单的SNR估算
    orig_snr = np.mean(original ** 2) / np.var(np.diff(original))
    proc_snr = np.mean(processed ** 2) / np.var(np.diff(processed))
    return 10 * np.log10(proc_snr / (orig_snr + 1e-10))


def calculate_peak_regularity(intervals: np.ndarray) -> float:
    """计算峰值规律性"""
    if len(intervals) == 0:
        return 0.0
    cv = np.std(intervals) / (np.mean(intervals) + 1e-10)
    return max(0.0, 1.0 - cv)


def analyze_heart_rate_trend(heart_rates: np.ndarray) -> Dict[str, Any]:
    """分析心率趋势"""
    if len(heart_rates) < 3:
        return {'trend': 'insufficient_data'}
    
    # 线性回归分析趋势
    x = np.arange(len(heart_rates))
    slope = np.polyfit(x, heart_rates, 1)[0]
    
    if abs(slope) < 0.1:
        trend = 'stable'
    elif slope > 0:
        trend = 'increasing'
    else:
        trend = 'decreasing'
    
    return {
        'trend': trend,
        'slope': float(slope),
        'trend_strength': abs(slope)
    }


def detect_heart_rate_outliers(heart_rates: np.ndarray) -> Dict[str, Any]:
    """检测心率异常值"""
    q1 = np.percentile(heart_rates, 25)
    q3 = np.percentile(heart_rates, 75)
    iqr = q3 - q1
    
    lower_bound = q1 - 1.5 * iqr
    upper_bound = q3 + 1.5 * iqr
    
    outliers = (heart_rates < lower_bound) | (heart_rates > upper_bound)
    
    return {
        'outlier_count': int(np.sum(outliers)),
        'outlier_ratio': float(np.sum(outliers) / len(heart_rates)),
        'outlier_indices': np.where(outliers)[0].tolist()
    }


def classify_heart_rate_stability(cv: float) -> str:
    """分类心率稳定性"""
    if cv < 0.05:
        return 'very_stable'
    elif cv < 0.1:
        return 'stable'
    elif cv < 0.2:
        return 'moderate'
    else:
        return 'unstable'


def generate_quality_recommendations(quality: QualityResults) -> List[str]:
    """生成质量改善建议"""
    recommendations = []
    
    if quality.quality_score < 0.5:
        recommendations.append("信号质量较差，建议检查传感器连接和信号采集环境")
    
    if quality.snr < 10:
        recommendations.append("信噪比较低，建议增强信号预处理或改善采集条件")
    
    if quality.peak_regularity < 0.7:
        recommendations.append("峰值检测不够规律，建议调整峰值检测参数")
    
    if quality.artifacts_ratio > 0.2:
        recommendations.append("信号伪影较多，建议加强伪影检测和去除")
    
    if not recommendations:
        recommendations.append("信号质量良好，可以进行进一步分析")
    
    return recommendations


def interpret_rmssd(rmssd: float) -> str:
    """解释RMSSD值"""
    if rmssd < 20:
        return "较低，可能表示副交感神经活动减少"
    elif rmssd < 50:
        return "正常范围"
    else:
        return "较高，表示良好的副交感神经活动"


def interpret_sdnn(sdnn: float) -> str:
    """解释SDNN值"""
    if sdnn < 50:
        return "较低，整体HRV较小"
    elif sdnn < 100:
        return "正常范围"
    else:
        return "较高，整体HRV较大"


def interpret_pnn50(pnn50: float) -> str:
    """解释pNN50值"""
    if pnn50 < 3:
        return "较低，副交感神经活动可能不足"
    elif pnn50 < 15:
        return "正常范围"
    else:
        return "较高，副交感神经活动活跃"


def interpret_lf_hf_ratio(ratio: float) -> str:
    """解释LF/HF比值"""
    if ratio < 1:
        return "副交感神经占主导"
    elif ratio < 2:
        return "自主神经平衡"
    else:
        return "交感神经占主导"


def assess_autonomic_nervous_system(hrv: HRVResults) -> Dict[str, str]:
    """评估自主神经系统状态"""
    assessment = {}
    
    # 基于多个指标的综合评估
    if hrv.rmssd > 30 and hrv.pnn50 > 10:
        assessment['parasympathetic'] = 'active'
    elif hrv.rmssd < 20 and hrv.pnn50 < 5:
        assessment['parasympathetic'] = 'low'
    else:
        assessment['parasympathetic'] = 'moderate'
    
    if hrv.lf_hf_ratio > 2:
        assessment['sympathetic'] = 'dominant'
    elif hrv.lf_hf_ratio < 1:
        assessment['sympathetic'] = 'low'
    else:
        assessment['sympathetic'] = 'balanced'
    
    # 整体评估
    if assessment['parasympathetic'] == 'active' and assessment['sympathetic'] == 'balanced':
        assessment['overall'] = 'healthy_balance'
    elif assessment['parasympathetic'] == 'low' and assessment['sympathetic'] == 'dominant':
        assessment['overall'] = 'stress_pattern'
    else:
        assessment['overall'] = 'moderate_function'
    
    return assessment


def print_analysis_summary(analysis: Dict[str, Any]):
    """打印分析摘要"""
    print("\n--- 分析摘要 ---")
    
    # 基本信息
    basic = analysis['basic_info']
    print(f"Pipeline: {basic['pipeline_name']}")
    print(f"信号长度: {basic['signal_length']} 点")
    print(f"检测峰值: {basic['peaks_detected']} 个")
    
    # 心率分析
    if 'heart_rate_analysis' in analysis and 'error' not in analysis['heart_rate_analysis']:
        hr = analysis['heart_rate_analysis']
        print(f"平均心率: {hr['mean_heart_rate']:.1f} BPM")
        print(f"心率稳定性: {hr['stability']}")
    
    # 质量分析
    if 'quality_analysis' in analysis and 'error' not in analysis['quality_analysis']:
        quality = analysis['quality_analysis']
        print(f"信号质量: {quality['quality_grade']} (评分: {quality['overall_score']:.3f})")
    
    # HRV分析
    if 'hrv_analysis' in analysis:
        hrv = analysis['hrv_analysis']
        print(f"HRV - RMSSD: {hrv['time_domain']['rmssd']:.1f} ms")
        print(f"自主神经评估: {hrv['ans_assessment']['overall']}")


if __name__ == "__main__":
    print("PPG结果分析模块")
    print("请使用 analyze_ppg_results() 和 visualize_ppg_results() 函数分析PPG处理结果")